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Abstract. Recently, with the advances in computational speed and availability there has 
been a growth in the number and resolution of fully 3-D hydrodynamical simulations. 
However, all of these simulations are purely hydrodynamical and there has been little at- 
tempt to include the effects of radiative transfer except in a purely phenomenological man- 
ner because the computational cost is too large even for modem supercomputers. While 
there has been an effort to develop 3-D Monte Carlo radiative transfer codes, most of these 
have been for static atmospheres or have employed the Sobolev approximation, which limits 
their applicability to studying purely geometric effects such as macroscopic mixing. Also 
the computational requirements of Monte Carlo methods are such that it is difficult to cou- 
ple with 3-D hydrodynamics. Here, we present an algorithm for calculating 1-D spherical 
radiative transfer in the presence of non-monotonic velocity fields in the co-moving frame. 
Non-monotonic velocity flows will occur in convective, and Raleigh-Taylor unstable flows, 
in flows with multiple shocks, and in pulsationally unstable stars such as Mira and Cepheids. 
This is a first step to developing fully 3-D radiative transfer than can be coupled with hydro- 
dynamics. We present the computational method and the results of some test calculations. 

1. Introduction 

The equation of radiative transfer (RTE) in spherical symmetry for moving media has been 
solved with a number of different metho ds, e.g. Monte Carl o calculation s (Magnan 197(i 
ICaroff et al.lll972t lAuer & Blerkomlll972l). Sobolev methods ( ICastoilll97Cil). the tangent rav 
metho d ('Mihalas. Kunasz. & Hu mmer"l976'). and the DOME method (Ikauschildt & Wehrs j 
Il99ll) . Today's state-of-the-art computer codes use iterative methods for the solution of the 
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RTE, based on the philosophy of operator splitting or operator perturbation JCannonlllQ?! 
IScharmei<fl98i . Following these ideas, different approximate A-operators for this "acceler- 
ated A-iteratio n " (ALI) metho d have been used successfully dOlson. Auer. & Buchleiin"987t 
lHamannl ll987>. meme?^l9^ and have been appli ed to the cons truction of non-LTE, ra- 
diative eq uilibrium models of ste llar atmospheres dWerneil fl987l) . Hauschildt (1992) and 
iHauschildt. Storzer. & Baro n (1994) have deve loped an opera tor splitting method based on the 
short-characteristic method (nOlson et al..l987;ii01son & Kunasz.1987.) to obtain the formal solu- 
tion of the special relativistic, spherically symmetric RTE along its characteristics and a band- 
diagonal approximation to the discretized A-operator. This method can be implemented very 
efficiently to obtain an accurate solution of the spherically symmetric RTE for continuum and 
line transfer problems using only modest amounts of computer resources. 

The main restriction of the co-moving frame (CMP) method discussed in lHauschildtl ( ll992l) IS 
a restriction to monotonic velocity fields. For monotonic velocity fields, the wavelength deriva- 
tive part of the CMF RTE can be posed as an initial value problem with the initial conditions 
set at small wavelengths for expanding media and at long wavelengths for contracting media. 
The initial value problem has to be solved by a fully implicit wavelength discretization (e.g., 
upwind schemes) in order to guarantee stability. In media with non-monotonic velocity fields, 
the wavelength derivative changes the structure of the equation so that it becomes a boundary 
value problem with boundary conditions at both short and long wavelengths. Since the equation 
is first order in wavelength at each spatial point there is only one boundary condition, whose 
wavelength sense depends on the local sign of the coefficient of the derivative. Therefore, if the 
wavelength-space is discretized to solve the CMF RTE, we need local upwind schemes in order 
to guarantee stability and to properly account for the presence of mixed boundary conditions in 
wavelength-space. 

In principle, it is possible to solve the RTE in the observer's frame, however, in that frame 
the emission and absorption processes are anisotropic and a detailed calculation is extremely 
complex. The main obstacle at high velocity, is that a prohibitively large number of wavelength 
and angle points are required to solve the observer's frame RTE. 

In this paper, we describe an operator splitting method to solve the spherical relativistic CMF 
RTE for arbitrary velocity fields. The method can be applied to a wide variety of astrophysical 
problems; e.g., atmospheres of pulsating stars, stellar atmospheres with shocks, multi-component 
novae, and supernova atmospheres. Although we present the method for the spherically symmet- 
ric 1-dimensional case, it can be extended to 3D geometry. Our approach calculates the formal 
solution along each characteristic independently (for known source function) in the combined 
radius-wavelength space. This has to be done in order to account for the wavelength and radial 
boundary conditions simultaneously during the formal solution stage of the operator splitting 
method. As the "approximate A operator" we use a block matrix that is tri-diagonal in wave- 
length space where each spatial (i.e., radial) block is itself a band matrix. This "wavelength 
tri-diagonal" approach yields superior convergence compared to a simpler "wavelength diago- 
nal" approach (which would correspond to a diagonal approximate A operator in static radiative 
transfer problems). The work presented here is a first step in an effijrt to develop methods for 
solving the CMF RTE in fully 3-D geometry, but even in its spherical version presented here will 
be useful for studying varying stars, such as Miras and Cepheids. 



2. Method 

In the following discussion we use notation of lHauschildtl lll992ll . Our starting point is the spher- 
ically symmetric form of the special relativistic, time independent (d/dt = 0) RTE, the restriction 
to plane parallel g eometry is straightforward. The calculation of the characteristics is identical to 
lHauschild3 ( 1 1 9921) and we thus assume that the characteristics are known. First, we will describe 
the process for the formal solution, then we will describe how we construct the approximate A 
operator. A* . 
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2.1. Formal solution 

We begin with Eq. 16 of lHauschildl ( ll992l) (see also Eq. 1 of lHauschildt & BaroDl2004l) : 

dAI , , .J 

-r + "-1-^ ^ ^li - V(i + ^ai)Ii (1) 
as OA 

where ds is a line element along a (curved) characteristic, Ii{s) is the specific intensity along the 
characteristic at point s >Q {s - Q denotes the beginning of the characteristic) and wavelength 
point A). The coefficient a/ is defined by 



r dr 



where p - v/c, j - ^\ - ffi and r is the radius. 77/ and xi the emission and extinction 
coefficients at wavelength /!/, respectively. 

Equation n describes the change of the intensity along a arbitrary characteristic though the 
medium. We define 

and discretize the wavelength derivative using a 3 point differencing formula (this can be made 
more general) to obtain: 

dh 

-r +ai \pu-\h-\ + piih + pu+\h+\\ = ??/ - Or/ + 4fl/)// (2) 
as 

where S\ - r]ilx\ is the source function at /I = Xi, dr - xds, and the pij are the discretization 
coefficients for the dAIi/dA derivatives (see Eqnsl^^Sjfor details). 

To obtain an expression for the formal solution, we rewrite Eq.|2]as 

dr 
with 

A XI ^ rji 
J I - — o / = — 

XI XI 
~ ai 

Si = -— [pi,i-ili-i + Pijli + Pu+ili+i\ 
XI 



(note that the index / denotes the wavelength point Ai) 

With this we obtain t he following expression for the formal solution (see also Eq. 14 in 
iHauschildt & Baronl2004 

//,/ = //- 1 ,/ exp(- At/_ 1 ) + 6Iij + 51; J (3) 

with the definitions 

SU,i - aijSi-xj + PijSij + YijSi+ij 

and 

Slij = aijSi-ij +j3ijSij 

The index i labels ffie (spatial) points along a characte ristic, the index / d enot es the wavelength 
point. The coefficients a, /?, and y, / are given in iHauschildJ ( 1 19921) and lOlson & Kunasj 
( 119871) . here they are calculated for a fixed wavelength for all points along a characteristic. S 
is a vector of known quantities (the old mean intensities and thermal sources). The S contain ffie 
effects of the velocity field on the formal solution and are given by 

c '^'■-1.' r r T r 1 

'J /-I,/ - [PlJ-lii-l,l-l + Pljii-lJ + Pl,l+lii-l,l+\\ 

Xi-lJ 

S i,i - --r- {pi.i-\h,i-\ + Pijh.i + Pi,i+\h.i+\\ 

Xi,l 
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This term provides the coupling of different wavelengths in the formal solution. Note that we 
have only used lir iear interpo l ation of the § terms and that this differencing scheme is dif- 
ferent than that of Hauschildt ( 1992). We have described this differencing scheme in detail in 
iHauschildt & Baron (■2004) . We have found that this differencing scheme is less diffusive in the 
co-moving frame and thus produces better line profiles for the case of small differential expan- 
sion. It gives identical results in the case of large global differential expansion (novae and super- 
novae) where the numerical diffusion in the co-moving frame is overwhelmed by the global line 
width. We have shown that this differencing scheme can become unstable under c ertain condi- 
tions, but can be stabilized in a straightforward manner (Hau schildt & Baronl2004i) . We present 
here the more complex ^ - 1.0, the generalization of the differencing scheme to the fully stable 
differencing scheme is obvious and can be found in Hauschildt & Baron (2004). 

If the velocity field is monotonically increasing or decreasing, then pij+i = or p/,/-! = for 
a stable upwind discretization of the wavelength derivative. In these cases, the problem becomes 
an initial value problem and can be solved for each wavelength once the results of the previous 
(smaller or longer) wavelength points are known. For non-monotonic velocity fields this is no 
longer the case and the formal solution needs to explicitly account for the wavelength couplings 
in both the blue and red directions. 

The formal solution is equivalent to the solution of one linear system for each characteristic. 
The rank of the system is «/ x n, where «/ is the number of wavelength points and «, is the number 
of points along the characteristic. For radial points, we have 3 < «, < 2«r - 1 points along 
each characteristic. The number of wavelength points n/ can be much larger, «/ x 1000 for the 
test cases presented later in this paper but «/ s; 300000 in full scale applications, so the rank of 
the systems can become large. Fortunately, the linear systems have a simple structure that allows 
us to use efficient methods for their numerical solution. 

For the construction of the system matrix for the formal solution along each characteristic it 
is useful to write the coefficients kij of the intensities as follows; 



o:i,iai-\APUi-\ 



Xi-U 

h-\,i ■ ki-ij = exp(-AT,_i) 



U,i : hi 



Xi-iJ 
/3ijai,iP 1,1-1 



Xi,i 

PijauPu 
Xi,i 

PijaijPu+i 
XiJ 

Jijaj+UPU-i 
Xi+iJ 

Jijaj+UPij+i 
Xi+iJ 



These expressions show the relatively simple matrix structure that can be exploited to solve for 
the mean intensities. 



2.1 .1 . Discretization of OAI/dA 



In order to ensure numerical stability, we use a local upwind scheme to discretize the wavelength 
derivative in the RTE: 
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- For 0; > 0: 



dAI 



All I - 

Ai - Ai-i 



- For ai < 0: 



dAI 
'dl 



All I - Ai+ili+i 



(4) 



(5) 



Ai - Ai+i 

Here, and in the following, we use a sorted wavelength grid with < Ai < Ai+i. The wave- 
length derivative is evaluated at a fixed spatial point along the characteristic. The coefficients p 
are then given by 



- For fl/ > 0: 



P 1,1-1 - - 



Pij = 



A,-i 



Ai - 
Ai 



Ai - 
Pl.M = 



- For fl/ < 0: 

Pu-i = 

Ai - Ai+i 
_ _ Ai+i 
Ai - Ai+i 

These coefficients depend on both the radial coordinate (because «/ is a function of r) and on 
the wavelengths (in the general case of a wavelength grid with variable resolution). Note that the 
direction of flow of information is determined by the sign of the coefficient a/ and not just on the 
velocity gradient. 



2.1.2. Boundary Conditions 

The spatial boundary conditions in the non-monotonic case remain the same as in the monotonic 
case, the incoming intensities at the spatial boundaries of each characteristic must be prescribed 
at every wavelength point. 

The wavelength boundary conditions are a bit more complicated, since now at every spatial 
point there is a wavelength boundary condition which must be determined by the local flow of 
information which is determined by the sign of each «/ along the characteristic. I.e., at each 
spatial point, one must determine whether information is flowing from blue-to-red or red-to-blue 
(in wavelength) and implement the proper boundary condition. 



2.1 .3. Structure of the system matrix 

We label the total n umber of wavelength p oints and the number of intersection points along 
a characteristic (see lHauschildt et aDll994 ) with n,. The total number of intensities that need to 
be determined is thus «, x n; per characteristic. Note that «, depends on the characteristic that is 
used in spherical symmetry and in the general 3-D case. There are two difi^erent ways to write 
the vector / of the specific intensities: 

1. "i-ordering": / - so that / is a vector of «/ vectors each of which has n, components. 

2. "1-ordering": / - so that / is a vector of n, vectors each of which has n; components. 

The position of the element (/, /) is therefore: 
1. "i-ordering": - (I - 1)«, + i and 
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2. "1-ordering": /;,■ — (i - l)n/ + I. 

In the following, we will denote these "block-vectors" with («,, «/) for both ordering schemes. 
If we write Eq.|3in matrix form, we obtain 

I ^AI + S + S 

where 5 is a (n,, «;) vector with the thermal and scattering source functions for each wavelength 
point «/ and radial point n, functions for each wavelength point, 5 is a (n,, «/) vector with the 
wavelength derivative information, and A is a (n, , «/) x (n,, «/) matrix. The row (/, I) of A has the 
following entries (the location of the element for i-ordering is also given): 



location 


index 


matrix element 


(i- 1,1-1) 


(/ - 2)n/ -H / - 1 




a -1,1 + 1) 


Irii + i — 1 




(i,l-l) 


(I - 2)n/ + i 


hi-i 


a -1,1) 


(I - l)n; -H / - 1 




(i,i) 


(I - l)n; + i 


hi 


(i,l+l) 


Irii + ; 


kij+i 


a + 1,1-1) 


(/ - 2)11; + i+l 


ki+ 1,1-1 


(i + 1,1 + 1) 


Irii + i + 1 


ki+ 1,1+1 



The total bandwidth of A in the i-ordering scheme is Zn, + i + 1 - ((I - 2)n,- + i — 1) - 2(n,- + 1). 
The 1-ordering scheme is symmetric to the i-ordering, thus the bandwidth in 1-ordering is simply 
2(«; -Hi). For large n;, i-ordering will require much smaller bandwidth, therefore, we consider 
only the i-ordering scheme here. 

To calculate the intensities along each characteristic we have to solve the hnear system 

{1-A)I^S+S (6) 

for each characteristic where / is the identity matrix. In the i-ordering scheme, the system matrix 
= / - A is thus a block-tridiagonal matrix with «/ blocks of n,- x n,- matrices. Fig. [J shows 
the general structure of the matrix y[. The / - 1 and / + 1 blocks are tridiagonal matrices, the 
/-block has the layout shown in Fig.|2 Note the sub-diagonal in the /-block is just - exp(-AT,_i). 
For the static case, the system degenerates to one with just the /-block non-zero allowing for 
direct recursive solution of the problem. For monotonic velocity fields either the / - 1 or the / + 1 
blocks are zero, again admitting recursive solution. In the case of non-monotonic velocity fields, 
all blocks may be non-zero and the system must be solved explicitly. It is convenient to call the 
tridiagonal matrix labeled / + 1 in Figure ^iM/^er, the lower diagonal matrix labeled /, diag, and 
the tridiagonal matrix labeled / - 1, sub. Then we can refer to e.g., the lower diagonal, diagonal, 
and upper diagonal of super as A "''^'"', B"''"^'', and C""'"''', respectively. 



2.1 .4. Solution of the linear systems 

The linear system of Eq. |6l has a relatively simple structure and can be s olved directly, e.g., 
by block tri-diagonal system solvers as described in lOolub & Loa n' ('1989). The problem with 
this approach is that the inverse of a tri-diagonal matrix is, in general, a full matrix. Therefore, 
the CPU time and memory requirements for direct solvers increase dramatically with increas- 
ing The special form of Eq.|5]and the sparseness of the blocks within the matrix A led us to 
investiga te iterative methods for the solution of Eq.|6] We examined the use of the "Rapido" al- 
gorithm jZurmiihl & Falkll986l page 194ff), however, as with all iterative linear system solvers, 
the eigenvalues of the matrix can (and are in many cases) be such that this method fails. For this 
work we also used standard band matrix solvers available in LAPACK and ESSE, which have the 
advantage that the pivots only have to b e calculated once per ALO iteration. Finally we used the 
general sparse solver package SuperLU dXiaove & D emme? '2003') which leads to large speedups 
over the other linear systems solvers we have tried (see below for a discussion of timing). 
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l+l 



Fig. 1. Structure of the system matrix. 



Fig. 2. Structure of the / block. 



2.2. Construction of A* 



The major difference between the monotonic velocity field problem and the non-monotonic ve- 
locity field is that we now have to construct a combined spatial-wavelength approximate A oper- 
ator for use in the operator splitting scheme. The basic equations of the operator splitting method 
remain unchanged, however In the discussion of the construction of A* it is useful to consider 
the "spatial" part of the A matrix and the "wavelength" part of the A matrix. The spatial part of 
the A matrix describes the transfer of photons in the radial coordinate and is, essentially, identical 
to the monotonic velocity field problem. The wavelength part of the A matrix describes the trans- 
port of photons in wavelength space due to the velocity field (or other non-coherent scattering 
processes). 

The A operator at wavelength point /, A;, has contributions from all wavelength points. This 
would lead to a matrix of order niXrir, where tir denotes the number of radial points. With a band- 
matrix form for the spatial A*'s, this would lead to a global band-matrix with significant storage 
requirements. Therefore, we derive a tri-diagonal approximation to the wavelength contributions 
to the global A* matrix. 

The construction of the spatial part of A* proceeds in exactly the same way as described in 

iHauschildt et al.. (.1994) . that is we assume a pulse of intensity of value unity is inserted into the 
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characteristic and by acting the A m atrix on the pulse we are able to construct and approximate 
lambda operator (ALO Ols on et al ][l987; Olson & Kunasz 1987). Even though we have made 
use of the fact that the formal solution can be written as a tridiagonal operator we shall show that 
each component of the ALO contains effects of both spatial and wavelength propagation of the 
global A matrix. 

We describe the construction of A* for arbitrary (spatial) bandwidth using the example of a 
tangential ray (core intersecting rays are a simple specialization of this case): The intersection 
points (including the point of tangency) are labeled from left to right, the direction in which the 
formal solution proceeds. For convenience, we label the ray tangent to shell / + 1 as /. Therefore, 
the ray / has 2/ + 1 points of intersection with discrete shells 1 . . . / + 1 . For each point k along 
the ray there is a "mirror point" A:,,, - 2i + \ - k. To compute row j of the discrete A-operator (or 
A-matrix), A,j, we sequentially label the intersection points of the ray / with the shell /' ("running 
index"), and define auxiliary quantities ^. The pulse is inserted at point ks, which is either k-\ 
or point k -Q 'm the case j — 1. It is convenient to define an X-factor as follows: 



Then 

^1/ = 



X 



1 _ 



^'kj-l ~ 



^kj-lkj 
^ ~ ^kj-l 
siih 



'^kj+l^k 



ks.l 



kj+l 

Then, propagating the pulse through the grid we obtain 
^k+l 

— ~ 1 

- k 

= + 1 



XiJ- 
Xjl 

X,i 

Then for Uj^ < we have 

W — l\ vdiag ^-Ir .sub W . a'''"S e' 1 

W _ ^, T,diag^-\r .diag i-i .supers .r/super^i -, 

W.'-i - ~ '^j.i-y LA;,/-!^ j-i.i-i "^jj-l i j-l.l ^ '^jj-l ^i,/J 
and for ap > we have 

t:' - n fjdiag.-lf.diag^i .supers -, 

With this we can construct A*'s with the full spatial bandwidth and tridiagonal in wavelength. 
However, in order to obtain good convergence, we also have to include wavelength dependent 
information in the global (spatial plus wavelength) A*. 

We can write the / component of the formal solution /fs in the form 

(/fs)/ = (A5)z * A,,,- 1 5,-1 + Ay5/ + A,,,+i5,+i (7) 
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Here, A;j_i, A/ /, and A/_/+i are the contributions to the A operator at wavelength point I origi- 
nating from the wavelength points / - 1, I, and I + 1. These matrices can be computed directly 
from the A'ji calculated above, simply by integrating over angle. Equation0shows explicitly the 
dependence of the mean intensities on the source functions of the neighboring wavelength points. 
Therefore, we can construct a global A* operator with the definition 

(A*S), = A*,,,_i5,-i + A*,jSi + A*i,mSm 

This leads to a block tri-diagonal global A*, where each spatial block is again a band matrix with 
the band- width given by the full spatial A*. This system can be solved efficiently either by direct 
solvers, using the same methods as discussed above for the formal solution. The convergence 
of the operator splitting method using the above A* is similar to that of static operator splitting 
methods with a (spatial) tri-diagonal ALO. Typically, 15-20 operator splitting iterations are re- 
quired to reach an accuracy of 10"^ for test cases with about 1000 wavelength points (see below). 
A wavelength-diagonal operator converges much more slowly, and an operator that ignores the 
wavelength dependence would require a maximum of «/ x niter iterations where niter is the number 
of iterations required to solve a monotonic velocity field problem for a single wavelength. 



3. Application examples 

As a first step we have implemented the method as a serial Fortran 95 program. This allows us 
to test the approach on problems with a relatively small number of wavelength points. Figure |3l 
describes the ste ps invo l ved in calcul ating the solution. Our basic test problem is similar to that 
discussed in lHauschildtU 19921) and in Ha uschildt & Baron (2004). We use a spherical shell with 
a grey continuum opacity parameterized by a power law in the continuum optical depth Tstd. The 
basic model parameters are 

1. Inner radius = lO'"' cm, outer radius rout = lO'^cm. 

2. Minimum optical depth in the continuum = 10""* and maximum optical depth in the 
continuum T™" - 10"^. 

std 

3. Grey temperature structure with Teff - 10 K. 

4. Outer boundary condition I^^ = and diffusion inner boundary condition for all wavelengths. 

5. Continuum extinction ;t'c = C/r^, with the constant C fixed by the radius and optical depth 
grids. 

6. Parameterized coherent & isotropic continuum scattering by defining 

Xc ^ ecKc + (l- ec)o-c 

with < Cc < I. Kc and cr^ are the continuum absorption and scattering coefficients. 

7. A parameterized spectral line with a rest wavelength of Aq - lOOOA and an intrinsic width 
of O.IA (equivalent to a width of 30km s"'), the line strength is parameterized by the ratio 
Xi(^o)/Xc, where is the line extinction coefficient. 

8. Parameterized line scattering defined analogous to the continuum scattering with a parameter 
e;. We assume complete redistribution for the line scattering. 



3.1. Tests for static and monotonic velocity fields 

This series of tests is designed to verify the correct operation of the code in the monotonic 
velocity case where we can compare it directly to our existing working code. 
The monotonic velocity field tests assume a linear velocity law of the form 

\f out 

where j3o = vq/c is the velocity at rout- 

Figurel^displays the flux transformed to the observer's frame for a linear velocity law defined 
above with vq - 1000 km s ' . The result is identical that produced by the purely monotonic code 
to the accuracy of our models. (We require that the scattering problem be solved to a relative 
accuracy of 10 '^.) We have also examined the moments in the co-moving frame and they too 
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For all wavelength 




For all characteristics 






build row 1 in 
the super-matrix 




\i 




write results to 

disk file as 
direct access 
records for bands 











Operator splitting iteration 



For all characteristics 



For all wavelength 




recall bands from 
disk file 






build 
supermatrix 






solve the linear 
system by 
any method 



I 



collect the 
specific intensities 
& build-up J's 

f 

Perform an OS/ALI 
step with the previously 
computed and saved ALO 



Fig. 3. Flowchart for the formal solution process. 

are identical to that obtained by the purely monotonic code. Thus, the formal solution and ALO 
solver and produce correct results in the regimes where we can test them directly. 

Table [2 shows the relative time for four different matrix solvers: 99% of the time is spent 
in solving the linear system in the formal solution, thus the effect of varying the matrix solvers 
is directly related to the code speed. Interestingly, SuperLU with the matrix stored in standard 
LAPACK band form is slower than LAPACK, where the matrix factors are saved to disk and re- 
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Matrix Solver 


Time 


LAPACK simple 


740s 


LAPACK witli recall 


468s 


SuperLU inefficient storage 


545s 


SuperLU efficient storage 


178s 



Table 1. The relative wall-clock time for a homology test with different matrix solvers. SuperLU 
with efficient storage leads significant speedup. 



called as needed; load and store is slower than I/O! In the recall case, along a given characteristic 
the ALO matrix is fixed and thus once it has been factored it need not be factored again, however 
there are too many factors to store in memory, since the problem size is larger than available 
memory, particularly when scaled to real problems requiring 300,000 wavelength points. Thus, 
the factors must be written to disk and read into the appropriate arrays as each characteristic is 
solved for. 

A large speedup comes from storing the SuperLU vector directly with no non-zero elements 
stored, even though this involves a considerable number of "if" statements in the construction 
loop. SuperLU can be called in a recall mode where the factors are stored similar to the LAPACK 
routines, but we expect to obtain much higher speedup by moving to SuperLUJDIST, the paral- 
lelized version of SuperLU and that wiU be discussed in a future paper 

3.2. Non-monotonic velocity fields 

In order to test our algorithm we have assumed the velocity structures shown in Figs.|5}|6l The 
first is a sine wave in zone number, such a structure could occur in a Mira or Cepheid variable star. 
The sine has been exponentially damped to mimic the effects of pho ton viscosity, whi ch should 
strongly damp the oscillations as the material becomes optically thin jMihalas & Mihalas.1984.) . 
The second, is just a illustrative set of piecewise continuous linear velocities, that might occur in 
e.g., the internal shock model of gamma-ray bursts. Figure0shows the continuum optical depth 
as a function of zone number for comparison. 

Figure |8] displays the observer's frame flux from the sine wave model for three choices of 
(ec, e/). The top panel is pure absorption in both the line and continuum, the middle panel has 
a pure absorptive line in a scattering continuum, and the bottom panel has a strongly scattering 
line in a scattering continuum. At first glance the line profile seems remarkably similar to the 
homologous, linear velocity law, except that there is a small dip near zero velocity. As scattering 
increases, the effect of the non-monotonic velocity law decreases. 

Figure|9]displays the observer's frame flux from the shock wave model with the same choices 
for (fc, e/) as in the sine wave model. The line profile is interesting, with only a small ripple in 
the flux profile at negative velocity. Surprisingly, this ripple does not seem to be washed out by 
scattering and hence is a feature of the velocity flow. 

In both cases (sine and shock) there are weak features that could be misidentified as real weak 
features if only homologous flows were considered. 



4. Conclusions 

We have presented a characteristics method of solving the full boundary value problem that 
occurs in both the wavelength and spatial dimension, using a full approximate lambda oper- 
ator in space, and tridiagonal in wavelength. The convergence of the operator is slow for the 
wavelength-diagonal operator, much improved for the wavelength-tridiagonal operator. Most of 
the computation time is spent in solving the linear system in t he formal solution. This c an be 
both sped-up and parallelized using e.g., the SuperLU package jXiaove & Demme]|l2003t) . The 
method we have presented is immediately useful for studying variable stars, and can be extended 
to 3-D radiation hydrodynamics and to studying other non-coherent scattering processes such 
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Fig. 4. The line profile in the observer's frame of the homologously expan ding (linear veloc- 
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ity law) model. The line is the results from the monotonia velocity code llHauschildt 
Hauschildt & Baronll2004 ) and the points are the results from the non-monotonic code for the 



case for =0.1 and ei = 10 . The results are identical. 

as partial redistribution, Compton scattering and neutrino transport in Raleigh-Taylor unstable 
flows that occur in core collapse supernovae. 
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Fig. 5. The velocity profile of the sin wave model. The exponential damping accounts for the 
damping efl'ects of photon viscosity. 
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Fig. 7. The continuum optical depth as a function of zone number for comparision. The optical 
depth scale is the same for both the sin wave model and the "shock" model. 
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